Criticality in polar fluids 
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■ A model of polar fluid is studied theoretically. The interaction potential, in addition to dipole- 

' dipole term, possesses a dispersion contribution of the van der Waals-London form. It is found 

that when the dispersion force is comparable to dipole-dipole interaction, the fluid separates into 
coexisting liquid and gas phases. The calculated critical parameters are in excellent agreement with 
Monte Carlo simulations. When the strength of dispersion attraction is bellow critical, no phase 
separation is found. 
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The phase behavior of simple molecular fluids was thought to be well understood since the early work of van der 
Waals over a century ago. It came, therefore, as a shock when the simulations in the early 90's could not locate-the 
anticipated liquid-gas coexistence in the simplest model of polar fluid, idealized by dipolar hard spheres Di/Stlcl. 
The shock was further augmented by the fact that our everyday experience with polar fluids, such as water, does 
not leave any doubt that they do posses both liquid and gas_phases. Furthermore, the theoretical methods such as, 
integral equationso and thermodynamic perturbation theorytm, developed over the hundred years since van der Waals 
first proposed his celebrated theory, all predicted such a phase separation. What could have gone wrong? 

The simulations of particles which interact by long ranged potentials, such as Coulomb or dipolar, are notorious for 
their difficulty. Thus, it is always possible that for some technical reason the simulations could not access the region 
i of instability. This, however, seems to be less and less likely in view of constantly improving computational power and 

algorithm design. The first hint of what was happening was provided by the original simulations of Weis, Levesque, 
. and Caillofl. These authors observed that at low temperatures, where the theories predicted phase separation into 
coexisting liquid and gas phases, the dipolar fluid became highly-structured, with particles forming chains of aligned 
dipoles. A further clue was obtained by van Leeuwen and Smitfl who simulated soft dipolar particles in which the 
softcore repulsion, 1/r 12 , was augmented by a variable isotropic attraction of van der Waals-London formLJ, 1/r 6 . 
^vq | By performing the Gibbs-ensemble Monte Carlo simulation van Leeuven and Smit came to the conclusion that a 
_^ ■ minimum strength of isotropic attraction was necessary in order to produce phase separation. The ubiquity of liquid- 
gas coexistence in real world can, therefore, be attributed to the presence of dispersion interactions between the 
molecules of polar fluids- A number of theories have been put forward to try to explain the unusual behavior found 
in computer simulationsutj, but none has, been able-tp completely resolve the mystery. Recently a new theory based 
on the pioneering ideas of Debye, HiickelJ, BjerrumE3, and OnsagerLj has been advanced to account for the absence 
of phase separation in DHSX-3. A similar approach has proven to be successful in the study of coulombic criticality in 
restricted primitive model (RPM) of electrolyte solutionsE3, predicting a coexistence curve in good agreement with 
experiments and computer simulationst 6 ! From the theoretical perspective, however, the system simulated by van 
Leeuven and Smito is significantly more complex than the RPM or the DHS. Besides the dipole-dipole interaction 
the theory must be able to deal with the softcore repulsion, as well as with the short ranged dispersion attraction. 
This model, therefore, provides a very stringent test of any theory of criticality in polar fluids. A successful theory 
• i— i , should account for both, the presence of phase separation for strong dispersion forces, as well as for its disappearance 
' when the isotropic coupling is reduced bellow the critical value. 

Our model, then, consists of N molecules of dipole moment p, inside a uniform medium of dielectric constant eo- 
- - > The interaction potential between any two particles is 

The dimensionless parameter A defines the strength of the isotropic interaction as it compares with dipole-dipole 
electrostatic potential. To proceed we associate with the model of "soft" molecules given by Eq. (flh an equivalent 
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system of "rigid" particles. The advantage of working with hard particles is that a number of analytic equations of 
state have been developed to study rigid molecules, in particular hard spheres. Fixing one molecule at the origin with 
the dipole moment aligned with the z-axis.-and a second particle at r with the dipolar orientation H, the potential 
in Eq. (|l]) can be separated in two partsEll. The repulsive term, U rp (r,Cl) — U(r;p 1 ,p 2 ) — U(r min ,fl min ) for 
\r\ < \r min \ and U rp (r,Sl) = for \r\ > \r min \; and the attractive term, U at {r,D) = U(r min ,n mm ) for \r\ < \r mm \ 
and U a t(r,Sl) = U(r;p l7 p 2 ) for |r| > |r mi „|. Here r min is the intermolecular vector and Sl min is the dipolar 
orientation for which the interaction potential is minimum. In the spirit of Weeks, Chandler, and Andersen (WCA)EJ 
we can associate with the soft particles given by Eq. (1) a system of rigid dipoles interacting by 



udd{r) = _ + J_ ( _ 3(P 1 -r-)(p 2 -r-) ) , 

jdd 



U aa {r) = 00 r < |a| , (2) 

where the distance of closest approach a((9; A) is determined by the condition, U rp (&, £l m in) = ksT . In view of 
the azimuthal symmetry, |a| = a(0; A) is only a function of the angle 9 between the dipole moment p 1 and the 
intermolecular vector r. To stress that this distance depends on the strength of the isotropic attraction, we have 
included A as an explicit parameter in a(0; A). The excluded volume region around the central dipole has the form of 
geoid, a flattened sphere, with the eccentricity a decreasing function of A (the highest eccentricity for A = 0). 

The reduced free energy density, / = (3F/V, of a fluid with the interaction potential Eq. || can be constructed 
as a sum of terms embodying the most relevant physical features, starting with the entropic ideal gas contribution, 
f ld = pln(pA 3 ) — p. Here p = N/V is the density of dipoles, — 1/fcsT, and A is the thermal wavelength. The mean 
energy of interaction between the dipoles can be calculated as 



jpvdW 



= \J dT X dT^)p{T 2 )U dd {v X ,T 2 ) , (3) 



where p(r) is the local density of particles and where, to simplify the notation, we have suppressed integration over 
the relative orientation of dipole moments. Since the distribution of particles and their orientation is uniform, only 
the first term of Eq. (|2|) contributes to the integral, with the second term averaging to zero. In order to simplify the 
integration we can replace the excluded volume region by a sphere with radius a s (A) = (a(0; A) + 2a(7r/2; A))/3. Since 
the eccentricity is quite small this is a good approximation. We have numerically checked that the error introduced 
by this approximation is no more than a few percent. Performing the integration, we find the reduced free energy due 
to dispersion interactions to be 

rvdW _ 27T A(/3*) 2 

1 3a3(A) T* ' 1 ' 

where we have defined the reduced (dimensionless) temperature and density to be, T* — kBTe^a 3 /p 2 and p* — pa 3 , 
respectively. 

Although on average the dipoles are distributed uniformly, with no particular preference for any specific orientation, 
the long ranged dipolar force produces strong correlations between the dipole moments of individual particles. This 
leads to a significant contribution to the total free energy. To obtain this correlational free energy, we fix one particle 
at the origin with its dipole moment aligned with the z-axis. The electrostatic potential that this dipole feels due to 
the presence of other molecules can be found from the solution of Laplace equation V 2 </> = 0. The boundary conditions 
require continuity of the potential, 4>i n (a s ) — 4>out{<is), and the displacement field, eo4>' in (a s ) = e(f>' out (a s ) , across the 
surface r = a s . The susceptibility of the outer dipoles to polarization by the electric field can be characterized by 
the repormalized dielectric constant e, the expression for which can be obtained from the Onsager's reaction field 
theoryEij, (e — eo)(eo + 2e)/e = Ait(3p 2 p. The Laplace equation can now be integrated to yield the electric field felt by 
the central dipole due to the other particles, 

E ° = J7T\o — 1 z ' 5 

e af(A) 2e + e 



The electrostatic free energy of the whole system is obtained from the Debye charging processE3L!j in which all the 
particles are charged simultaneously from ze ro ,, to their final dipolar strength. The integration can be done explicitly 
yielding the correlational free energy densitycJl23, 
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with ip{ u ) = e i u )/ e o = \ (1 + u ) + |v9 + 2m + u 2 , and m = Amp* jT*. Finally, the contribution to the total free 
energy due to hardcore repulsion can be approximated by the Carnahan-Starling (CS) forma, f hc = pg(rf), where 
g{rj) = 77(4 — 2>rf)/(l — rf) 2 and r\ = 7raf(A)p/6 is the volume fraction occupied by dipoles. 

Combining all of the contributions, the total free energy of dipolar fluid becomes / = f ld + f hc + f vdw + f dd . It 
is a simple matter to see that for any A this free energy violates the convexity requirement when the temperature is 
lowered below the critical value T C (A). The critical parameters are plotted in Fig. 1. It is evident that the agreement 
with the simulations is quite good. For 0.3 < A < 1, the critical densities are close to the ones obtained in the 
simulations. 
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FIG. 1. Critical density and critical temperature as a function of A: the dashed curve is for the theory without association; 
the two solid curves a and b are, for the limiting cases h(n) ~ 1.52165... and h(n) ~ 1.19055..., respectively. Solid circles are 
the Monte Carlo data from Ref[fl]. Inset shows the average cluster size as a function of A. 



For A = 1 the critical temperature is about 30% too high. However, for smaller isotropic couplings, A < 0.6, this 
discrepancy almost completely disappears. It is not surprising that for A ~ 1 the agreement is the poorest. After 
all, the theory presented above is intrinsically mean-field, and as such neglects the thermal fluctuations which are the 
strongest in the vicinity of the critical point. For short ranged potentials, the effect of fluctuations is to significantly 



depress the critica 
T c by 15% to 20% 



temperature. Even for such relatively simple system as "argon" mean-field theories overestimate 
itB. We might expect that as the relative strength of isotropic short ranged interaction decreases, 
compared to the long ranged dipolar force, the effect of fluctuations on the critical temperature will also diminish. 
This, indeed, is the case, see Fig. 1. We note that for Awl the phase separation is the result of the competition 
between the attractive and the excluded volume, CS, interactions. For small A's, on the other hand, the instability 
is purely electrostatic and is only weakly affected by the hardcore contribution to the total free energy. 

For A < 0.3 the simulations fail to find any coexistence, while our theory predicts phase separation all the way 
down to A = 0. What is responsible for this discrepancy? In view of the previous discussion, we do not have to look 
too far for an explanation. It was observed in simulations that for A < 0.3, the low temperature thermodynamics was 
dominated by polymer-like chains of aligned dipoles. On the other hand, the theory presented above is intrinsically 
linear, and does not take into account the strong non-linear effects associated with the formation of clusters. To 
account for these, while preserving the theory's linear structure, we shall explicitly allow for chains composed of 
n associated dipoles. Based on the Monte Carlo simulations we shall ignore other possible geometries of clusters. 
This, actually, is a non-trivial approximation, since at zero temperature the energy favors compact clustersli-l The 
simulations, however, find that for finite temperatures linear chains predominate. Thus, as a first approximation we 
shall only consider polymer-like clusters. 

The density of chains containing n monomers will be denoted by p n , with the density of free, unassociated, dipoles 
p\. The particle conservation leads to p = Y^=i n Pn: where the total density of dipoles is p = N/V. This is the 
so called "chemical", as opposed to "physical" picture of statistical mechanics of strongly interacting particles. It is 
important to keep in mind that no new approximations are being introduced into the theory, but only a change of 
perspective. Clearly if, as it happens at high temperatures, the formation of clusters is unlikely, the theory should 
inform us of this fact by predicting low cluster density. As-was noted by Onsager "what we remove from one page of the 
ledger would be entered elsewhere with the same effect"ll3. In fact it is possible to directly map this chemical picture 
onto exact thermodynamic perturbation theoryEl The new perspective is especially useful at low temperatures, where 
the perturbation theory is very difficult to construct, but where the "chemical" view point provides a particularly 
clear perspective of the physical phenomena. We can now proceed to make the first approximation. In the spirit of 
Bjerrun£3, we shall assume that the clusters interact only weakly with each other and with the free dipoles. Thus, as 
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a second approximation, we shall treat clusters as ideal non-interacting species with all of the interactions confined 
to free dipoles. This approximation has proven to work quite well in the theory of electrolytes, in particular allowing 
for a precise location of the critical pointta. 

The entropic free energy due to thermal motion of dipoles and clusters can be approximated by the ideal gas 
contribution, f ld = J2n°=i Pn m (y°nA 3 ™/Cn) — Pn- In the limit of low temperatures, the internal partition function of 
an n cluster, £„, can be evaluated for chainlike configurations to yieldOoO 

( 7r 3/2 2 21/18 T .5/2^»-l f(n-lWn)1 



where, 



h{n) 



n [^)(n) - ^/,( 2 ) (1)] + 2[^ l \n) - ^ (1)]1 4/3 

(8) 
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The ip^ (n) and (n) are the polygamma functions of the first and second order, respectively. For free dipoles (i = 1 . 
The hardcore contribution can, once again, be approximated by the Carnahan-Starling form, / = Y^n=i PnQiv)^ 
where r\ = ira^.(\)p/6 is the total volume fraction occupied by dipoles and chains. The mean-field and the electrostatic 
contributions to the free energy are given by Eqns. (Q) and (||), with the total density p replaced by the density of 
free dipoles p\. The distribution of cluster densities can be obtained through minimization of the total free energy 
/ = / + f hc + f vdw + f dd , subject to constraint of particle conservation. This reduces to the law of mass action, 
p n = npi, relating the chemical potential of free dipoles p\ — df/dpi, to the chemical potential of n clusters 
Pn = df /dp n . Substituting the free energy, the low of mass action yields 

Pn = Q n pn e Pn^+{n-l)9{v) ; „ > 2 ; (9) 

where the excess chemical potential is p ex = df dd /dp\. In spite of its simple appearance the law of mass action is 
actually an infinite set of coupled algebraic equations which determine the distribution of chain densities {p n }- To solve 
these equations is a difficult numerical task. In order to proceed we observe that h(n) is a uniformly increasing function 
bounded by the two limiting values, h(2) = 1.19055... and h(oo) = 1.52165...; i.e. 1.19055... < h(n) < 1.52165. ..Vn. 
Use of these bounds allows us to delimitate the critical parameters. Approximating h(n) by one of the limiting values, 
Eq. (H) can be summed explicitly, X^^LzJ 1 / 3 " = P~ Pit reducing to a single algebraic equation for p\. Substituting 
the root of this equation back into Eq. (g) leads to the distribution of cluster sizes. For large A's, separation in two 
coexisting phases is once again found. In fact, the critical densities and temperatures are identical to our earlier 
calculation done in the absence of cluster formation. As A decreases, the critical temperature also diminishes, and 
dipolar association becomes relevant. We find that below certain critical value, 0.07 < A c < 0.35, the phase separation 
disappears, see Fig. 1. Van Leeuwen and Smit could not locate phase separation for A < 0.3, the value which is very 
close to our upper bound. We note, however, that the lower bound is probably more realistic since, as can be seen 
from the inset to Fig 1, the average chain size is still quite small in the region of criticality. 

We conclude that the theory presented above captures the essential physics of criticality in polar fluids. A number 
of issues still need to be addressed. As it stands, the theory does not take into account interactions between the 
chains. How these residual forces can be introduced into the theory is far from obvious. We must stress, however, 
that the naive approach of alkjwing each one of the monomers inside the chains to interact with unrenormalized 
potential, Eq. ([!]), is incorrecto. This would double count the interactions, which are partially taken care of by the 
dipolar association, and would lead to phase separation for all values of A. 

The nature of criticality in polar fluids is also of great interest. Although there can be little doubt that the 
Stockmayer fluid (A = 1) belongs to the Ising universality class, in view of our results, this no longer seems to be so 
evident for dipolar fluids with A « 0.3. The fact that the mean-field theory can account so well for the location of 
the critical point, suggests that for A < 0.5 the fluctuations play only a marginal role. One might, therefore, expect 
to see mean-field critical exponents, or a crossover from the mean-field to the Ising universality in a very narrow 
neighborhood of the critical point. As is the case for electrolytes, the criticality in dipolar systems should remain a 
challenge for the foreseeable future. _ 

While this paper was going through the refereeing process, new simulationsE3 have once again raised a possibility of 
the phase separation of DHS at sufficiently low temperatures and densities. Unlike, the usual gas-liquid transition, 
it has been argued, however, that this novel phase separation is driven by the proliferation of Y-like bifurcations^! 
The two coexisting phases are found to have respectively a large and a small concentration of Y-like defectsEl 
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